Cell type specific markers in WT Arabidopsis roots¶
Working directory should contain root atlas Seurat object, which is available on GEO under record GSE152766:¶
Root_Atlas.rds (or get it by running through notebook 2, 3, 4, 5, 6, 7, 8, 9 & 10)
Working directory must contain subdirectories,supp_data¶
supp_data/ should contain the files, which are available on Github (github.com/Hsu-Che-Wei/COPILOT):¶
features.tsv.gz
Curated_Markers.csvIn [1]:
rm(list=ls())
# Set the working directory to where folders named after the samples are located.
# The folder contains spliced.mtx, unspliced.mtx, barcodes and gene id files, and json files produced by scKB that documents the sequencing stats.
setwd("/scratch/AG_Ohler/CheWei/scKB")
In [1]:
library(tidyverse)
library(Seurat)
library(RColorBrewer)
library(future)
#for 200gb ram
options(future.globals.maxSize = 200000 * 1024^2)
In [3]:
# Load atlas
rc.integrated <- readRDS("./Root_Atlas.rds")
In [4]:
# Simple QC label
rc.integrated$celltype.anno <- gsub("Putative Quiescent Center", "Quiescent Center", rc.integrated$celltype.anno, ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
order <- c("Quiescent Center", "Stem Cell Niche", "Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Pericycle", "Phloem", "Xylem", "Procambium", "Unknown")
palette <- c("#9400d3", "#DCD0FF", "#5AB953", "#BFEF45", "#008080", "#21B6A8", "#82B6FF", "#0000FF","#FF9900","#E6194B", "#9A6324", "#FFE119","#EEEEEE")
rc.integrated$celltype.anno <- factor(rc.integrated$celltype.anno, levels = order[sort(match(unique(rc.integrated$celltype.anno),order))])
color <- palette[sort(match(unique(rc.integrated$celltype.anno),order))]
In [5]:
# Plot atlas with celltype annotation
options(repr.plot.width=8, repr.plot.height=6)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", cols=color)
In [6]:
# Plot atlas with time zone annotation
options(repr.plot.width=8, repr.plot.height=6)
time_plt <- DimPlot(rc.integrated,
group.by = "time.anno",
order = c("Maturation","Elongation","Meristem"),
cols = c("#DCEDC8", "#42B3D5", "#1A237E"))
time_plt
FindMarkers for each cell type¶
In [7]:
# Set the identity of interest to cell type
Idents(rc.integrated) <- "celltype.anno"
# Set the expression matrix used for DE to batch-corrected and scaled one
DefaultAssay(rc.integrated) <- "integrated"
In [8]:
# run differentially expression (DE) analysis using ROC method, select randomly 10000 cells for each cell type
Clust_Markers <- FindAllMarkers(rc.integrated,
logfc.threshold=log(2),
min.diff.pct = 0.25,
max.cells.per.ident = 10000,
only.pos=T,
test.use="roc")
In [9]:
# Load in gene ids and names
feature_names <- read_tsv("./supp_data/features.tsv.gz", col_names = c("gene", "Name", "Type")) %>%
select(-Type) %>%
distinct()
In [10]:
# Merge the DE output with corresponding gene id and name
Clust_Markers <- left_join(Clust_Markers, feature_names)
In [11]:
# For each cell type, select the top 3 genes according to AUC score
Clust_Markers %>% group_by(cluster) %>% top_n(3, myAUC)
In [12]:
Clust_Markers$cluster_gene <- paste(Clust_Markers$cluster, Clust_Markers$gene, sep="_")
In [13]:
# store all results
All_clust_markers <- Clust_Markers
# subset based on AUC 0.75
Clust_Markers <- filter(Clust_Markers, myAUC>=0.75)
In [14]:
# Calculate the frequency of genes being detected as marker
(times_per_gene <- Clust_Markers %>%
ungroup() %>%
group_by(gene) %>%
tally())
In [15]:
Clust_Markers <- left_join(Clust_Markers, times_per_gene)
In [16]:
summary(Clust_Markers$n)
In [17]:
# Load libraries for plotting
library(cowplot)
library(ComplexHeatmap)
library(circlize)
library(GeneOverlap)
library(gprofiler2)
library(ggrepel)
library(ggplot2)
In [18]:
options(repr.plot.width = 16, repr.plot.height = 14)
markers_sel <- select(Clust_Markers, gene, cluster)
markers_list <- split(markers_sel, f=markers_sel$cluster)
#this makes list from long df of gene lists - TARGET is what we want to keep
markers_list <- lapply(markers_list, function(x) x[names(x)=="gene"])
# convert each sublist into character and eliminate duplicates
markers_list <- lapply(markers_list, function(x) as.character(unique(x$gene)))
## GeneOverlap
# number of integrated features
genome_size <- 17513L
#compare all lists
gom.self <- newGOM(markers_list, markers_list, genome.size=genome_size)
int <- getNestedList(gom.self, "intersection")
int_matrix <- getMatrix(gom.self, "intersection")
p.val <- getMatrix(gom.self, "pval")
JC <- getMatrix(gom.self, "Jaccard")
# log of p.val for intersection
p.val_log <- -log10(p.val + 1e-200)
olap <- Heatmap(p.val_log,
name = "-log10_pval",
col = colorRamp2(c(0, 200),
c("beige", "red")),
column_title = "Number of shared marker genes",
cluster_rows = T,
cluster_columns = T,
use_raster= FALSE,
show_column_names = TRUE,
show_row_names = TRUE,
show_row_dend = TRUE,
clustering_distance_rows = "pearson",
clustering_distance_columns = "pearson",
show_column_dend = TRUE, cell_fun = function(j, i, x, y, width, height, fill) {grid.text(sprintf("%.0f", int_matrix[i, j]), x, y, gp = gpar(fontsize = 10))
})
# padding - bottom, left, top, right
draw(olap, padding = unit(c(15, 5, 5, 10), "mm"), heatmap_legend_side = "left")
In [19]:
# Load GO terms for Arabidopsis thaliana
cluster_GO <- gost(markers_list, organism = "athaliana", correction_method = "fdr", significant = F, multi_query = F)
cluster_GO_df <- cluster_GO[[1]]
cluster_GO_sig <- filter(cluster_GO_df, p_value<=0.01)
# top terms for each cluster
cluster_GO_sig %>%
filter(source=="GO:BP", intersection_size>=4) %>%
group_by(query) %>%
top_n(5, wt = -p_value) %>%
arrange(desc(p_value)) -> top_GO
GO_n <- cluster_GO_sig %>%
filter(source=="GO:BP", intersection_size>=4) %>%
group_by(term_id) %>%
tally() %>%
arrange(desc(n))
GO_n <- dplyr::rename(GO_n, "n_clusters"=n)
cluster_GO_sig_n <- left_join(cluster_GO_sig, GO_n)
# get all terms for the top ones so that all clusters have values
top_GO_all <- filter(cluster_GO_df, term_id %in% top_GO$term_id)
#spread and plot
top_GO_sel <- select(top_GO_all, query, p_value, term_id, term_name)
spread_GO <- spread(top_GO_sel, key = query, p_value)
spread_GO[is.na(spread_GO)] <- 1
spread_GO_m <- as.matrix(-log10(spread_GO[3:ncol(spread_GO)]))
rownames(spread_GO_m) <- spread_GO$term_name
In [20]:
options(repr.plot.width = 18, repr.plot.height = 20)
GO_hm <- Heatmap(spread_GO_m,
name = "-log10_pval",
heatmap_legend_param = list(title_position="topcenter", color_bar = "continuous"),
col = colorRamp2(c(0, 10),
c("beige", "#e31a1c")),
cluster_rows = T,
cluster_columns = T,
use_raster= FALSE,
show_column_names = TRUE,
show_row_names = TRUE,
show_row_dend = TRUE,
show_column_dend = TRUE,
clustering_distance_rows = "pearson",
clustering_distance_columns = "pearson",
row_names_gp = gpar(fontsize = 12))
# padding - bottom, left, top, right
draw(GO_hm, padding = unit(c(15, 15, 5, 80), "mm"), heatmap_legend_side = "left")
In [21]:
# Load the known markers
(used.markers <- read_csv("./supp_data/Curated_Markers.csv"))
In [22]:
used.markers$cluster_gene <- paste(used.markers$Celltype, used.markers$Locus, sep="_")
In [23]:
# Combine DE output to see whether the known markers are captured
Clust_Markers <- mutate(Clust_Markers, pct.diff=pct.1-pct.2)
Clust_Markers <- arrange(Clust_Markers, desc(pct.diff)) %>%
group_by(cluster) %>%
mutate(pct.diff_rank=dplyr::row_number()) %>%
arrange(desc(avg_diff)) %>%
mutate(avg_diff_rank=dplyr::row_number()) %>%
arrange(desc(myAUC)) %>%
mutate(myAUC_rank=dplyr::row_number()) %>%
mutate(combined_rank_raw=(pct.diff_rank + avg_diff_rank + myAUC_rank)/3) %>%
arrange(combined_rank_raw) %>%
mutate(combined_rank=dplyr::row_number()) %>%
select(-combined_rank_raw) %>%
mutate(known_marker=cluster_gene %in% used.markers$cluster_gene) %>%
arrange(combined_rank)
Clust_Markers
In [24]:
# Store the data
write.csv(Clust_Markers, "./supp_data/Atlas_celltype_markers_ROC.csv", row.names=F)
In [25]:
# Plotting the expression pattern of known markers
used.markers.to.plt <- used.markers
used.markers.to.plt$cluster <- used.markers.to.plt$Celltype
used.markers.to.plt$cluster <- factor(used.markers.to.plt$cluster, levels= c("Quiescent Center", "Stem Cell Niche", "Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Pericycle", "Phloem", "Xylem", "Procambium"))
used.markers.to.plt <- arrange(used.markers.to.plt, cluster)
used.markers.to.plt
ct_to_plt_Names <- rev(used.markers.to.plt$Gene)
ct_to_plt_genes <- rev(used.markers.to.plt$Locus)
breaks <- rev(ct_to_plt_genes)
labels <- rev(ct_to_plt_Names)
In [26]:
rc.integrated$celltype.anno <- fct_rev(rc.integrated$celltype.anno)
Idents(rc.integrated) <- "celltype.anno"
In [27]:
options(repr.plot.width=14, repr.plot.height=8)
dot <- DotPlot(object = rc.integrated, features = ct_to_plt_genes, cols = c("white", "red")) + RotatedAxis() + ylab("") +
scale_x_discrete(breaks=breaks, labels=labels)
dot
In [28]:
# top genes to plot
novel_markers <- Clust_Markers %>%
group_by(cluster) %>%
top_n(2, -combined_rank)
novel_markers$cluster <- factor(novel_markers$cluster, levels= c("Quiescent Center", "Stem Cell Niche", "Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Pericycle", "Phloem", "Xylem", "Procambium"))
novel_markers <- arrange(novel_markers, cluster)
novel_markers
novel_markers_sel <- select(novel_markers, cluster, gene, Name)
ct_to_plt_Names <- rev(novel_markers_sel$Name)
ct_to_plt_genes <- rev(novel_markers_sel$gene)
In [29]:
(breaks <- rev(ct_to_plt_genes))
In [30]:
(labels <- rev(ct_to_plt_Names))
In [31]:
options(repr.plot.width=10, repr.plot.height=6)
dot <- DotPlot(object = rc.integrated, features = ct_to_plt_genes, cols = c("white", "red")) + RotatedAxis() + ylab("") +
scale_x_discrete(breaks=breaks, labels=labels)
dot
Specific genes by cell type and developmental stage¶
In [32]:
table(rc.integrated$time.celltype.anno)
In [33]:
### combination of time and cell type anno
rc.integrated$cell_time_final <- rc.integrated$time.celltype.anno
In [34]:
rc.integrated$cell_time_final <- gsub("Meristem", "Meristematic", rc.integrated$cell_time_final, ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
rc.integrated$cell_time_final <- gsub("Elongation", "Elongating", rc.integrated$cell_time_final, ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
rc.integrated$cell_time_final <- gsub("Maturation", "Mature", rc.integrated$cell_time_final, ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
# Meristematic-Putative Quiescent Center
rc.integrated$cell_time_final <- gsub("Meristematic-Putative Quiescent Center", "Meristematic-Quiescent Center", rc.integrated$cell_time_final, ignore.case = FALSE, perl = FALSE,
fixed = T, useBytes = FALSE)
In [35]:
table(rc.integrated$cell_time_final)
In [36]:
Idents(rc.integrated) <- "cell_time_final"
DefaultAssay(rc.integrated) <- "integrated"
cell_time_Clust_Markers <- FindAllMarkers(rc.integrated,
logfc.threshold=log(2),
min.diff.pct = 0.25,
max.cells.per.ident = 10000,
only.pos=T,
test.use="roc")
In [37]:
cell_time_Clust_Markers
In [38]:
cell_time_Clust_Markers %>%
group_by(cluster) %>%
tally()
In [39]:
cell_time_Clust_Markers %>%
filter(myAUC>=0.75) %>%
group_by(cluster) %>%
tally()
In [40]:
# store all
all_cell_time_Clust_Markers <- cell_time_Clust_Markers
In [41]:
# keep only markers 0.75 and above
cell_time_Clust_Markers <- filter(cell_time_Clust_Markers, myAUC>=0.75)
In [42]:
(times_per_gene <- cell_time_Clust_Markers %>%
ungroup() %>%
group_by(gene) %>%
tally())
In [43]:
cell_time_Clust_Markers <- left_join(cell_time_Clust_Markers, times_per_gene)
In [44]:
cell_time_Clust_Markers <- left_join(cell_time_Clust_Markers, feature_names)
In [45]:
cell_time_Clust_Markers$cluster_gene <- paste(cell_time_Clust_Markers$cluster, Clust_Markers$gene, sep="_")
Rank AUC dataset for celltype/dev combo¶
In [46]:
cell_time_Clust_Markers <- mutate(cell_time_Clust_Markers, pct.diff=pct.1-pct.2)
cell_time_Clust_Markers
In [47]:
cell_time_Clust_Markers <- arrange(cell_time_Clust_Markers, desc(pct.diff)) %>%
group_by(cluster) %>%
mutate(pct.diff_rank=dplyr::row_number()) %>%
arrange(desc(avg_diff)) %>%
mutate(avg_diff_rank=dplyr::row_number()) %>%
arrange(desc(myAUC)) %>%
mutate(myAUC_rank=dplyr::row_number()) %>%
mutate(combined_rank_raw=(pct.diff_rank + avg_diff_rank + myAUC_rank)/3) %>%
arrange(combined_rank_raw) %>%
mutate(combined_rank=dplyr::row_number()) %>%
select(-combined_rank_raw) %>%
arrange(combined_rank)
cell_time_Clust_Markers
In [48]:
options(repr.plot.width = 16, repr.plot.height = 14)
markers_sel <- filter(cell_time_Clust_Markers, myAUC>=0.75) %>% select(gene, cluster)
markers_list <- split(markers_sel, f=markers_sel$cluster)
#this makes list from long df of gene lists - TARGET is what we want to keep
markers_list <- lapply(markers_list, function(x) x[names(x)=="gene"])
# convert each sublist into character and eliminate duplicates
markers_list <- lapply(markers_list, function(x) as.character(unique(x$gene)))
## GeneOverlap
# number of integrated features
genome_size <- 17513L
#compare all lists
gom.self <- newGOM(markers_list, markers_list, genome.size=genome_size)
int <- getNestedList(gom.self, "intersection")
int_matrix <- getMatrix(gom.self, "intersection")
p.val <- getMatrix(gom.self, "pval")
JC <- getMatrix(gom.self, "Jaccard")
# log of p.val for intersection
p.val_log <- -log10(p.val + 1e-200)
olap <- Heatmap(p.val_log,
name = "-log10_pval",
col = colorRamp2(c(0, 200),
c("beige", "red")),
column_title = "Number of shared marker genes - AUC 0.75+",
cluster_rows = T,
cluster_columns = T,
use_raster= FALSE,
show_column_names = TRUE,
show_row_names = TRUE,
show_row_dend = TRUE,
clustering_distance_rows = "pearson",
clustering_distance_columns = "pearson",
show_column_dend = TRUE, cell_fun = function(j, i, x, y, width, height, fill) {grid.text(sprintf("%.0f", int_matrix[i, j]), x, y, gp = gpar(fontsize = 10))
})
# padding - bottom, left, top, right
draw(olap, padding = unit(c(15, 5, 5, 10), "mm"), heatmap_legend_side = "left")
In [49]:
markers_sel <- filter(cell_time_Clust_Markers, myAUC>=0.75) %>% select(gene, cluster)
markers_list <- split(markers_sel, f=markers_sel$cluster)
#this makes list from long df of gene lists - TARGET is what we want to keep
markers_list <- lapply(markers_list, function(x) x[names(x)=="gene"])
# convert each sublist into character and eliminate duplicates
markers_list <- lapply(markers_list, function(x) as.character(unique(x$gene)))
cluster_GO <- gost(markers_list, organism = "athaliana", correction_method = "fdr", significant = F, multi_query = F)
cluster_GO_df <- cluster_GO[[1]]
cluster_GO_sig <- filter(cluster_GO_df, p_value<=0.01)
# top terms for each cluster
cluster_GO_sig %>%
filter(source=="GO:BP", intersection_size>=4) %>%
group_by(query) %>%
top_n(5, wt = -p_value) %>%
arrange(desc(p_value)) -> top_GO
GO_n <- cluster_GO_sig %>%
filter(source=="GO:BP", intersection_size>=4) %>%
group_by(term_id) %>%
tally() %>%
arrange(desc(n))
GO_n <- dplyr::rename(GO_n, "n_clusters"=n)
cluster_GO_sig_n <- left_join(cluster_GO_sig, GO_n)
# get all terms for the top ones so that all clusters have values
top_GO_all <- filter(cluster_GO_df, term_id %in% top_GO$term_id)
#spread and plot
top_GO_sel <- select(top_GO_all, query, p_value, term_id, term_name)
spread_GO <- spread(top_GO_sel, key = query, p_value)
spread_GO[is.na(spread_GO)] <- 1
spread_GO_m <- as.matrix(-log10(spread_GO[3:ncol(spread_GO)]))
rownames(spread_GO_m) <- spread_GO$term_name
In [50]:
options(repr.plot.width = 18, repr.plot.height = 20)
GO_hm <- Heatmap(spread_GO_m,
name = "-log10_pval",
heatmap_legend_param = list(title_position="topcenter", color_bar = "continuous"),
col = colorRamp2(c(0, 10),
c("beige", "#e31a1c")),
cluster_rows = T,
cluster_columns = T,
use_raster= FALSE,
show_column_names = TRUE,
show_row_names = TRUE,
show_row_dend = TRUE,
show_column_dend = TRUE,
clustering_distance_rows = "pearson",
clustering_distance_columns = "pearson",
row_names_gp = gpar(fontsize = 12))
# padding - bottom, left, top, right
draw(GO_hm, padding = unit(c(15, 15, 5, 80), "mm"), heatmap_legend_side = "left")
In [51]:
cell_time_Clust_Markers$clus_to_split <- cell_time_Clust_Markers$cluster
cell_time_Clust_Markers <- separate(cell_time_Clust_Markers, clus_to_split, into=c("timezone", "celltype"), sep="-")
In [52]:
cell_time_Clust_Markers$cluster <- as.character(cell_time_Clust_Markers$cluster)
cell_time_Clust_Markers$timezone <- factor(cell_time_Clust_Markers$timezone, levels=c("Mature", "Elongating", "Meristematic"))
ct_levels <- names(table(rc.integrated$celltype.anno))
cell_time_Clust_Markers$celltype <- factor(cell_time_Clust_Markers$celltype, levels=ct_levels)
cell_time_Clust_Markers <- arrange(cell_time_Clust_Markers, timezone) %>% arrange(celltype)
cell_time_Clust_Markers
In [53]:
write.csv(cell_time_Clust_Markers, "./supp_data/Atlas_cell_time_combined_ROC_Clust_Markers_20200527.csv", row.names=F)
Selected Markers for Cell type and developmental stage specific dot plot¶
SCR, MYB36 and CASP1 were selected for the endodermis because they are highlighted later in the manuscript and are highly enrichened. Note that MYB36 is expressed in all developmental stages of endodermis, but is highest in the elongation zone.
For others, I took the top 20 combined rank and gave preference for those that were specific (i.e. n<=1) or where known markers
In [54]:
(Markers_to_plot <- read_csv("./supp_data/Atlas_cell_time_combined_ROC_Clust_Markers_20200527_to_plot.csv") %>% filter(plot=="y"))
In [55]:
(c_levels <- unique(Markers_to_plot$cluster))
In [56]:
ct_to_plt_Names <- Markers_to_plot$Name
ct_to_plt_genes <- Markers_to_plot$gene
In [57]:
table(rc.integrated$cell_time_final)
In [58]:
rc.integrated$cell_time_final <- factor(rc.integrated$cell_time_final, levels=rev(c_levels))
Idents(rc.integrated) <- "cell_time_final"
table(rc.integrated$cell_time_final)
In [59]:
(breaks <- rev(ct_to_plt_genes))
In [60]:
(labels <- rev(ct_to_plt_Names))
In [61]:
options(repr.plot.width=14, repr.plot.height=10)
DotPlot(object = rc.integrated, features = ct_to_plt_genes, cols = c("white", "red")) + RotatedAxis() +
scale_x_discrete(breaks=breaks, labels=labels) + xlab("")
In [62]:
options(repr.plot.width=14, repr.plot.height=10)
dot <- DotPlot(object = rc.integrated, features = ct_to_plt_genes, cols = c("white", "red")) + RotatedAxis() + ylab("") + theme(axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank()) +
scale_x_discrete(breaks=breaks, labels=labels) + xlab("")
dot
In [63]:
options(repr.plot.width=14, repr.plot.height=10)
dot <- dot + geom_rect(xmin = 0.5, xmax = 1.5, ymin = 0.5, ymax = 1.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 1.5, xmax = 2.5, ymin = 1.5, ymax = 2.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 2.5, xmax = 5.5, ymin = 2.5, ymax = 5.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 5.5, xmax = 8.5, ymin = 5.5, ymax = 8.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 8.5, xmax = 11.5, ymin = 8.5, ymax = 11.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 11.5, xmax = 14.5, ymin = 11.5, ymax = 14.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 14.5, xmax = 17.5, ymin = 14.5, ymax = 17.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 17.5, xmax = 20.5, ymin = 17.5, ymax = 20.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 20.5, xmax = 23.5, ymin = 20.5, ymax = 23.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 23.5, xmax = 26.5, ymin = 23.5, ymax = 26.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 26.5, xmax = 29.5, ymin = 26.5, ymax = 29.5, alpha = 0,fill=alpha("grey",0), color="black") +
geom_rect(xmin = 29.5, xmax = 32.5, ymin = 29.5, ymax = 32.5, alpha = 0,fill=alpha("grey",0), color="black")
dot
In [64]:
ann_df <- tibble(timezone=as.character(Markers_to_plot$timezone),
celltype= as.character(Markers_to_plot$celltype)) %>%
mutate(row=dplyr::row_number(),
col=rep("A")) %>%
group_by(celltype) %>%
mutate(ct_row=dplyr::row_number()) %>%
ungroup() %>%
mutate(Label=celltype)
ann_df$Label[ann_df$ct_row!=2] <- " "
ann_df$Label[ann_df$celltype=="Stem Cell Niche"] <- "Stem Cell Niche"
ann_df$Label[ann_df$celltype=="Quiescent Center"] <- "Quiescent Center"
ann_df$row <- factor(ann_df$row)
ann_df$timezone <- factor(ann_df$timezone, levels=c("Meristematic", "Elongating", "Mature"))
In [65]:
ann_df
In [66]:
options(repr.plot.width=1, repr.plot.height=10)
(time_ann_p <- ann_df %>% ggplot(aes(x=col, y=rev(row), fill=timezone)) + geom_tile() + scale_fill_manual(values = c("#DCEDC8", "#42B3D5", "#1A237E")) + theme_nothing())
In [67]:
options(repr.plot.width=2, repr.plot.height=10)
ann_df %>% ggplot(aes(x=col, y=rev(row), fill=timezone)) + geom_tile() + scale_fill_manual(values = c("#DCEDC8", "#42B3D5", "#1A237E")) + theme_minimal()
In [68]:
options(repr.plot.width=1, repr.plot.height=10)
ann_df$celltype <- factor(ann_df$celltype, levels = order[sort(match(unique(ann_df$celltype),order))])
color <- palette[sort(match(unique(ann_df$celltype),order))]
(cell_ann_p <- ann_df %>% ggplot(aes(x=col, y=rev(row), fill=celltype)) + geom_tile() + scale_fill_manual(values = color) + theme_nothing())
In [69]:
options(repr.plot.width=1, repr.plot.height=10)
cell_ann_p | time_ann_p
In [70]:
options(repr.plot.width=2.2, repr.plot.height=10)
(text_ann_p <- ann_df %>% ggplot(aes(x=col, y=rev(row), label=Label)) + geom_text(size = 5) + theme(plot.margin = margin(0, 0, 35, 35, "cm")) + theme_nothing())
In [71]:
options(repr.plot.width=4.5, repr.plot.height=10)
(c_ann <- text_ann_p | (cell_ann_p | time_ann_p))
In [72]:
options(repr.plot.width=14, repr.plot.height=10)
combined_marker_plot <- text_ann_p + cell_ann_p + time_ann_p + dot + patchwork::plot_layout(widths = c(0.8,0.2,0.2,4), nrow=1)
combined_marker_plot
In [73]:
options(repr.plot.width=9, repr.plot.height=6)
order <- c("Quiescent Center", "Stem Cell Niche", "Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Pericycle", "Phloem", "Xylem", "Procambium", "Unknown")
order <- rev(order)
palette <- c("#BD53FF", "#DED3FE", "#5AB953", "#BFEF45", "#008080", "#21B6A8", "#82B6FF", "#0000FF","#FF9900","#E6194B", "#9A6324", "#FFE119","#EEEEEE")
palette <- rev(palette)
rc.integrated$celltype.anno <- factor(rc.integrated$celltype.anno, levels = order[sort(match(unique(rc.integrated$celltype.anno),order))])
color <- palette[sort(match(unique(rc.integrated$celltype.anno),order))]
cell_dim <- DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", cols=color)
cell_dim[[1]]$layers[[1]]$aes_params$alpha = 1 # can change w/ this if you like. I decided not to so side anno would match
cell_dim
In [74]:
options(repr.plot.width=9, repr.plot.height=6)
time_plt <- DimPlot(rc.integrated,
group.by = "time.anno",
order = c("Meristem","Elongation","Maturation"),
cols = c("#1A237E", "#42B3D5", "#DCEDC8"))
time_plt
In [75]:
options(repr.plot.width=16, repr.plot.height=18)
(cell_dim + time_plt) / combined_marker_plot + patchwork::plot_layout(nrow=2, heights = c(1, 1.4), guides = 'collect')
In [76]:
ggsave("./supp_data/Atlas_cell_dev_marker_dot.pdf", width=16, height=18)
In [77]:
sessionInfo()